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ABSTRACT 


One  of  the  more  significant  error  sources  in  using  the  Global  Positioning  System 
for  attitude  determination  of  aircraft  and  earth  orbiting  satellites  is  multipath  interference. 
The  GPS  signal  is  reflected  fi'om  the  spacecraft  structures  resulting  in  a  composite  signal 
that  differs  firom  the  direct  line-of-sight  signal  in  terms  of  phase  and  amplitude.  The 
resulting  phase  errors  can  have  a  significant  impact  on  the  attitude  determination  error. 
This  thesis  will  describe  the  GPS  signal  generation  and  provide  tools  to  simulate  the 
signal.  A  Kalman  filter  was  developed  to  track  the  phase  changes  in  the  simulated  GPS 
signals  and  its  performance  described  in  the  absence  and  presence  of  multipath. 
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1. 


INTOODUCTION 


A.  BACKGROUND 

In  order  for  a  spacecraft  to  perform  its  mission,  it  must  be  able  to  orient  itself  in 
relation  to  either  the  sun,  the  earth  or  both  within  an  inertial  reference  frame.  The  primary 
system  used  to  accomplish  this  is  the  spacecraft  attitude  determination  and  control 
subsystem  (ADCS).  The  ADCS  uses  external  references  to  determine  the  vehiele’s 
absolute  attitude  with  respect  to  some  fked  inertial  reference  frame.  Traditionally,  these 
external  references  have  included  the  Sun,  the  Earth’s  infrared  (IR)  horizon,  the  local 
magnetic  field  direction,  and  the  stars.  Recent  technological  developments  have  added  a 
new  tool  to  the  external  reference  arsenal  -  the  Global  Positioning  System  (GPS). 

Typical  traditional  ADCS  sensor  characteristics  are  shown  below  in  Table  1 . 
(Larson  and  Wertz) 


SENSOR 

TYPICAL  PERFORMANCE 

RANGE 

WEIGHT  (KG) 

POWER  (W) 

Inertial  measurement  unit 

(gyros  &  accelerometers) 

Gyro  drift  rate  =  0.003°/hr  to 

r/hr 

3  to  25 

10  to  200 

Sun  Sensors 

Accurac5F=0.005°  to  3° 

0.5  to  2 

0to3 

Star  Sensors  (scarmers  & 

mappers) 

0.0003°  to  0.01° 

3  to  7 

5  to  20 

Horizon  Sensors 

O 

O 

o 

O 

2  to  5 

5  to  10 

Magnetometer 

0.5°  to  3° 

0.6  to  1.2 

<  1 

Table  1  Typical  ADCS  Sensors 


In  comparison,  the  Trimble  TANS  Vector  attitude  determination  system  using 
GPS  is  capable  of  an  accuracy  of  0. 1°  with  a  mass  of  2.2  kg  (GPS  receiver  and  antennas) 
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and  a  power  requirement  of  7.5  W.  Additionally,  the  Trimble  system  also  provides  three- 
dimensional  orbital  position  information;  which  other  sensors  are  incapable  of  providing. 

As  described  by  Cohen  (1996),  attitude  determination  using  GPS  is  achieved  by 
determining  the  relative  positions  of  multiple  antennas,  mounted  to  the  vehicle,  based  on 
differences  in  the  carrier  phase  measurements  at  the  different  antennas.  The  basic  premise 
is  that  the  GPS  satelhte  is  so  distant,  relative  to  the  receiving  antenna  separations,  that 
arriving  wavefronts  can  be  considered  as  effectively  planar  as  shown  below  in  Figure  1. 


Figure  1  Attitude  Geometry 

The  phase  of  a  signal,  ^ ,  is  a  function  of  the  wavelength,  X,  and  distance  traveled, 
d: 

d 

(f>  =  ln:— 

The  wavelength  is  a  function  of  the  speed  of  light,  c,  and  the  frequency,  f  The 
distance  traveled  is  a  function  of  the  speed  of  light  and  the  time  difference  of  arrival 
between  antennas.  At: 
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d  =  cM 

Therefore  a  relationship  between,  time  and  phase  can  be  shown  to  be 
cAt 

(p  -  Ittj  =  Ik-j-j  -  Ar 

The  incoming  signal  will  arrive  at  the  antenna  closest  to  the  transmitter  before 
reaching  the  most  distant  antenna.  By  measuring  the  difference  in  carrier  phase  between 
the  antennas,  a  receiver  can  determine  the  relative  range  between  any  two  anteimas.  The 
measured  differential  phase,  A(p,  for  basehne  /  and  satellite  /,  is  proportional  to  the  dot 
product  of  the  baseline  vector,  3c ,  and  the  line  of  sight  unit  vector  to  the  transmitter,  s . 

As  shown  in  Figure  2,  the  GPS  receiver  measures  only  the  fractional  part  of  the 
differential  phase.  The  integer  component,  k,  must  be  resolved  by  independent  means. 

The  differential  phase  can  therefore  be  expressed  as: 

A<Py  =  s^ -x-k^j+Vy 

where  is  the  noise  in  the  measurement.  Given  the  differential  phase 

measurements  between  the  different  antennas,  the  spacecraft  attitude  can  then  be  resolved 
as  shown  by  Cohen  (1995).  Filtering  this  measurement  will  be  the  focus  of  this  thesis. 

The  dominant  noise  source  in  differential  phase  measurements  is  multipath,  which 
occurs  when  the  signal  arriving  at  the  antenna  consists  of  the  reflected  signals  in  addition 
to  the  line-of-sight  source.  The  reflected  signal  is  phase  shifted  with  respect  to  the  original 
transmission  and  appears  as  additive  noise  at  the  anteima.  Because  the  antenna  locations 
are  different,  the  multipath  signature  at  each  antetma  is  unique  and  will  not  be  common  to 
every  antenna.  Multipath  accounts  for  more  than  90%  of  the  total  error  budget  in  carrier 
phase  measurements.  (Lightsey) 
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B.  PURPOSE 


This  thesis  will  discuss  the  generation  of  the  GPS  carrier  signal  which  is  used  in 
making  differential  phase  measurements.  Models  will  be  developed  for  simulating  the 
multipath  environment  and  the  GPS  signal.  A  Kalman  filter  which  tracks  the  phase  in  the 
signal  for  use  in  the  differential  phase  measurements  will  be  described.  Finally,  results  of 
the  Kalman  filter’s  use  in  a  multipath  environment  will  be  presented. 

Chapter  n  will  discuss  the  GPS  Signal  Generation  and  present  a  model  of  the 
signal.  The  multipath  interference  will  be  examined  and  modeled  in  Chapter  IE.  Chapter 
IV  provides  the  derivation  of  the  Kalman  filter  for  tracking  the  phase  of  a  GPS  signal. 
Chapter  V  will  show  the  results  when  the  filter  is  used  in  tracking  the  GPS  phase.  The 
final  chapter  will  summarize  the  findings  and  present  other  ideas  for  minimizing  multipath 
and  suggestions  for  future  investigation. 
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n.  GPS  SIGNAL  GENERATION 


The  GPS  signal  is  generated  as  shown  in  Figure  3.  The  GPS  satellite  transmits 
two  carrier  frequencies,  LI  (primary  frequency)  and  L2  (secondary  frequency).  These 
carrier  frequencies  are  modulated  by  spread  spectrum  codes  with  a  pseudorandom  noise 
(PRN)  sequence  unique  to  each  satellite.  The  LI  frequency  is  modulated  by  the 
navigation  message  along  with  two  PRN  codes:  the  coarse/acquisition  (C/A)  code  and  the 
precision  (P)  code.  The  GPS  L2  signal  is  modulated  by  one  of  the  codes  using  Binary 
Phase  Shift  Keying  (BPSK).  The  GPS  LI  signal  is  modulated  using  both  codes  by  an 
unbalanced  quadrature  phase  shift  key  (QPSK)  modulator  which  is  actually  two  BPSK 
modulators  with  different  amplitudes.  This  development  closely  follows  that  given  by 
Kaplan  (1996). 

A.  PN  CODE  GENERATION 

The  PRN  codes  used  in  the  generation  of  the  GPS  signal  are  non-maximum  length 
codes  also  known  as  Gold  codes.  These  codes  are  generated  from  an  n-bit  shift  register 
generator.  This  register  has  an  initial  state  or  fill.  After  each  clock  cycle  the  output  of  a 
specified  cell  is  used  as  the  input  to  a  modulo-2  adder  (exclusive-or).  The  output  of  the 
modulo-2  adder  is  then  fed  back  into  a  specified  stage  of  the  register.  As  an  example, 
consider  the  shift  register  shown  in  Figure  4.  This  register  uses  the  3^“*  and  5*  stages  as 
the  input  to  the  modulo-2  adder  and  feeds  the  result  back  into  the  first  stage.  The 
polynomial  used  to  describe  the  design  of  linear  code  generators  is  of  the  form 

1+Z^'' 

where  x‘  means  the  output  of  the  i*  cell  of  the  shift  register  is  used  as  the  input  of 
the  modulo-2  adder  and  the  1  means  the  output  of  the  adder  is  feed  in  to  stage  1  The 
polynomial  for  the  register  in  Figure  4  is  then  1  +  +  X^  If  this  register  was  initially 

filled  with  all  I’s  at  time  0,  the  resulting  register  for  the  first  ten  clock  cycles  would  be  as 
shown  in  Table  2.  The  binary  signal  output  is  taken  from  the  last  register  and  in  this  case 
would  be  11111000110. 
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Information 


Figure  3  GPS  LI  and  L2  Signal  Generation 


Clock 


Output 


Shift  Register 

Figure  4  PRN  Code  Shift  Register  for  Example  Polynomial  1 


time 
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4 

5 
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1 

1 

0 

0 

0 

1 

1 

0 

1 
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1 
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Register 
2  3  4  5 

1111 
1111 
0  111 
0  0  11 
0  0  0  1 
10  0  0 
110  0 
0  110 
10  11 
110  1 
1110 


Table  2  PRN  Code  Generation  Example  Results  for  Polynomial  1 


The  Matlab  m-fiinction  gen  poly.m  can  be  used  to  generate  the  PRN  sequence  for 
any  user  defined  polynomial  and  register. 


B.  C/A  CODE  GENERATION 

The  GPS  coarse/acquisition  (C/A)-code  is  a  Gold  code  generated  by  two  10-bit 
shift  registers,  G1  and  G2.  The  polynomial  for  each  is  given  below  (Kaplan,  1996) 

G1  =  1+XHx^° 

G2  =  1  +  +  X®  +  X*  + 

Both  C/A-code  registers  have  an  initial  fill  of  all  I’s.  Each  code  is  derived  and  the 
result  of  G2  is  then  delayed  by  the  satellite  PRN  code  number.  This  delay  effect  is 
obtained  by  the  exclusive-or  of  the  selected  positions  of  the  G2  register.  The  C/A-code  is 
then  obtained  by  the  exclusive-or  of  the  G1  and  delayed  G2  signals.  An  example  for 
satellite  vehicle  1  is  shown  in  Figure  5. 

The  code  tap  positions  and  initial  code  sequences  are  given  in  Table  3.  The  m-files 
gl.m  and  g2.m  generate  the  codes  for  a  specified  satellite  vehicle  number. 
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Figure  5  C/A-Code  Generator  for  Satellite  Vehicle  1  with  C/A  Code  Tap  Selection  2Q6 

(from  Kaplan) 
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14 

1772 

4343 

15 

809 

257 

15 

1775 

4343 

16 

9010 

258 

16 

1776 

4343 

17 

104 

469 

17 

1156 

4343 

18 

205 

470 

18 

1467 

4343 

19 

306 

471 

19 

1633 

4343 

20 

407 

472 

20 

1715 

4343 

21 

508 

473 

21 

1746 

4343 

22 

609 

474 

22 

1764 

4343 

23 

103 

509 

23 

1063 

4343 

24 

406 

512 

24 

1706 

4343 

25 

507 

513 

25 

1743 

4343 

26 

608 

514 

26 

1761 

4343 

27 

709 

515 

27 

1770 

4343 

28 

8010 

516 

28 

1774 

4343 

29 

106 

859 

29 

1127 

4343 

30 

307 

860 

30 

1453 

4343 

31 

308 

861 

31 

1625 

4343 

32 

409 

862 

32 

1712 

4343 

Table  3  GPS  Code  Phase  Assignments  (from  GPS-ICD-2000) 
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test  CA.m  is  a  script  file  to  test  C/A  code  generation.  This  program  calculates 
the  first  10  C/A-chips  and  converts  the  result  to  octal  numbers  that  match  the  First  10 
C/A-Chips  (Octal)  of  the  GPS-ICD-2000  table  shown  in  Table  3 

For  example  the  first  10  chips  of  the  PRN  5  C/A-code  is  1001011011  which  is 
1133  in  octal. 

C.  P-CODE  GENERATION 


The  GPS  Precision  (P)  Code  is  generated  by  PRN  sequences  using  four  12  bit  shift 
registers  whose  polynomials  and  initial  states  are  shown  below: 


Register 

Polynomial 

Initial  State 

XIA 

1  +  X^  +  X^+  x^+  x^+  x^+  x^^+  x^^  +  x^^ 

001001001000 

XIB 

1  +  x^  +  x^+  x^  +  x^+  x^  +  x’  +  x**+  x^+  x‘”+  x“+  x‘^ 

010101010100 

X2A 

1  +  x^  +  x^+  X‘+  x^+  x’+  x*+  x^+  x^‘’  +  x^^  +  x^^ 

100100100101 

X2B 

1  +X2  +  X^  +  X^  +  X*  +  X®  +  X*^ 

010101010100 

Table  4  P-Code  Generator  Polynomials  and  Initial  States 


The  P-code  is  7  days  in  length  at  a  chipping  rate  of  10.23  MBps.  The  7  day 
sequence  is  the  Modulo-2  sum  of  two  subsequences  XI  and  X2i;  their  lengths  are 
15,345,000  chips  and  15,345,037  chips  respectively. 

XI  itself  is  generated  by  the  Modulo-2  sum  of  the  output  of  two  12-stage  registers 
(XI A  and  XIB)  short  cycled  to  4092  and  4093  chips  respectively.  When  the  XI A  short 
cycles  are  counted  to  3750  the  XI  epoch  is  generated.  This  occurs  every  1.5  seconds, 
after  15,345,000  chips  of  the  XI  pattern  have  been  generated.  The  XI  period  is  defined 
as  3750  XI A  cycles  (15,345,000  chips)  which  is  not  an  integer  number  of  XIB  cycles.  To 
accommodate  this  the  XIB  shift  register  is  held  in  the  final  state  (chip  4093)  of  its  3749* 
cycle  for  343  chips.  The  X2  sequence  is  generated  in  a  similar  manner,  however  the  X2 
period  is  defined  to  be  37  chips  longer  than  the  XI  period  in  order  to  force  the  X2A  and 
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X2B  epochs  to  precess  with  respect  to  the  Xla  and  XIB  epochs.  This  precession 
continues  until  the  end  of  the  7  day  period. 

The  Y-code*  is  used  when  the  anti-spoofing  mode  of  operation  is  activated.  It  is 
generated  in  the  same  fashion  as  the  P-code,  however  the  polynomial  generating  equations 
are  classified.  (lCD-GPS-2000) 

xla.m,  xlb.m,  x2a.m  and  x2b.m  are  used  to  calculate  the  PN  codes.  The  m-file 
testj^.m  is  used  to  test  P-code  generation. 

D.  DATA  GENERATOR 

The  output  of  the  data  generator  is  the  navigation  message,  Z)(r) .  The  message 
includes  satellite  vehicle  ephemerides,  system  time,  clock  bias  information,  status 
messages,  and  C/A  to  P-code  handover  information.  The  navigation  message  is  generated 
at  50  bps  and  is  added  to  the  P  and  C/A-codes  to  modulate  both  the  LI  and  L2  signals. 
(ICD-GPS-200) 

E.  BPSK  MODULATOR 

After  the  various  codes  are  generated  they  are  then  used  to  modulate  the  carrier 
signal  using  binary  phase  shift  keying,  BPSK,  resulting  in  the  LI  or  L2  signal.  The 
Simulink  program  used  to  simulate  the  BPSK  modulator  is  shown  below  in  Figure  6. 


*  The  Y-code  and  P-code  are  sometimes  discussed  together  and  are  referred  to  as  P(Y)-code  meaning  P  (or 
Y)  code. 
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Outport 

Figure  6  BP SK  Simulink  Diagram 

The  PRN  code  is  input  via  the  inport.  It  is  sent  through  a  zero-order-hold  and 
relay  function  to  generate  a  non-return  to  zero  (NRZ)  signal,  which  is  multiplied  by  the 
sine  wave,  resulting  in  a  binary  phase  shifted  signal. 
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BPSK  Simulation 


0  2  4  6  8  10 


Figure  7  BPSK  Simulation  Results 

Figure  7  shows  an  example  of  the  BPSK  simulation  results.  This  figure  shows  the 
results:  NRZ  signal  in  the  top  graph,  the  sine  wave  in  the  second  graph,  the  product  of  the 
two,  or  the  BPSK  signal  in  the  third  graph,  and  a  composite  of  the  data  and  resulting 
signal  in  the  bottom  graph. 

F.  LI  SIGNAL 

The  P-code  at,  10.23  x  10^  chips  per  second,  and  C/A-code,  at  1.023  x  10^  chips 
per  second,  are  combined  with  the  50  bits  per  second  data  and  are  then  used  to  modulate 
the  L 1  fi'equency  ( 1 54  •  /<,  =  1 575.42  MHz )  as  shown  in  Figure  3 .  This  results  in  an 
Unbalanced  Quadri-Phase  Shift  Keying  (UQPSK)  modulation  with  the  data  bits  added  to 
the  ranging  codes,  the  C/A-code  signal  lagging  the  P-code  signal  by  90°;  and  the  C/A- 
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code  signal  power  nominally  exceeding  the  P-code  signal  power  by  3  dB.  (Struza) 
Mathematically,  this  can  be  represented  by: 

L^{co^t)  =  A\P(t)  0  Z)(/)]cos(fi?,?)  +  V2^,[G(/f)  0  i>(/)]  sin(<»/) 
where  Pit)  is  the  P-code,  Dit)  is  the  data,  and  G(r)  is  the  C/A-code. 

G.  L2  SIGNAL 

The  L2  jfrequency  ( 120-  /,  =  1227.6  Aff/z)  is  modulated  by  either  the  exclusive-or 
of  the  P-code  and  data,  the  P-code,  or  the  exclusive-or  of  the  C/A  code  and  data,  as 
selected  by  the  control  segment.  Mathematically,  this  can  be  represented  by: 

L^ico^t)  -  ■O(O]cos(02O 

L^ico^t)  =  4[P(0]cos(o20 
Liicojt)  =  4[G(0]cos(fi>20 


H.  SUMMARY 

The  GPS  signal  consists  of  two  separate  signals,  LI  and  L2.  These  signals  are 
generated  using  BPSK  and  unbalanced  QPSK  modulation.  The  modulation  signal  is 
generated  using  the  P  and  C/A  codes  to  which  the  data  signal  is  added.  The  Matlab  fQes 
initm,  sig  gen.m  and  my  glb.m  are  used  to  initialize  and  generate  a  sample  of  the  codes 
used  to  modulate  the  signal.  These  codes  are  then  used  by  the  Simulink  programs 
LI  sim.m  and  L2  sim.m  to  generate  a  sample  of  the  signal.  In  the  next  chapter  the 
multipath  environment  will  be  described. 
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in.  MULTIPATH 


A.  INTRODUCTION 

This  chapter  will  discuss  the  multipath  environment  that  impacts  the  use  of  GPS 
for  attitude  determination  and  will  present  a  method  to  model  multipath  interference  using 
Matlab. 

Multipath  occurs  when  the  signal  arrives  at  the  antenna  from  reflected  surfaces  and 
from  the  line  of  sight  sources.  A  simple  case  with  a  single  reflected  path  signal  is  shown  in 
Figure  8. 


The  reflected  signal  is  phase  shifted  with  respect  to  the  original  transmission  and 
appears  as  additive  noise  at  the  antenna.  When  antenna  locations  are  different,  the 
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multipath  signature  at  each  antenna  is  unique  and  the  error  in  the  phase  measurement  is 
not  common.  (Lightsey)  Multipath  interference,  due  to  environmental  factors,  will  impact 
the  signal,  however  this  effect  will  be  common  to  each  received  signal. 

For  the  developed  model  it  will  be  assumed  that  the  multipath  environment  will  be 
static.  In  reality  the  multipath  environment  will  be  var5nng  due  to  the  satellites  attitude 
and  position  relative  to  the  GPS  satellites  changing  over  time. 

B.  MATHEMATICAL  MODEL 

In  the  absence  of  the  multipath,  the  input  signal  to  a  GPS  receiver  for  the  LI 
QPSK  signal  is  given  by; 

5(0  =  cos(eo  t  +a{t)f) 

where  Ac  is  the  received  signal  amplitude  and  a(t)  =  0,1,2  or  3  is  the  pseudo¬ 
random  code  waveform  that  phase  modulates  the  carrier.  (Kumar)  In  the  presence  of  N 
multipaths,  in  addition  to  the  direct  line  of  sight  path,  the  input  signal  is  given  by 

Sn,(t)  =  a^s(t)  +  a,s(t  -  r, )+. .  .+ajv‘S(^  -  ) 

where  and  r,.  denote  the  amplitude  and  delay  of  the  i**'  multipath  for 
i=l,2,. . .  ,N.  Graphically,  this  can  be  represented  as  shown  in  Figure  9. 
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Typically  for  BPSK  and  QPSK  signals  the  carrier  frequency  operates  at  a  much 
higher  frequency  than  the  modulation  rate.  As  a  result,  a  small  value  of  delay,  t  ,  will 
result  in  large  phase  differences.  When  the  multipath  signal  with  a  large  phase  difference 
is  added  to  the  direct  path  there  may  be  destructive  adding  which  can  significantly  impact 
the  received  signal, 

C.  SIMULATION  MODEL 

In  the  model  mpath.m,  three  vectors  are  used  to  develop  the  signal 

delay  -1  by  vector  of  delayed  samples.  The  first  entry  is  the  undelayed 
signal,  the  next  entry  is  the  signal  delayed  by  one  sample  period.  The  length  of  this  vector 
is  the  maximum  value  in  the  tau  vector  and  is  initialized  as  all  zeros. 

tau  -  vector  of  delay  times.  Each  entry  is  1  minus  the  delay  time  in  sample 
periods,  dt,  and  is  used  to  index  the  delay  vector.  The  first  entry  is  1  and  represents  a 
delay  of  zero. 

alpha  -  vector  of  amplitudes.  Each  entry  is  the  amplitude  of  the 
corresponding  delayed  signal.  The  first  entry  is  a  1  in  order  to  return  the  undelayed  line 
of  sight  signal. 

As  an  example,  the  BPSK  signal  shown  in  Figure  10,  with  N=3  multipath  signals, 
will  be  used  to  demonstrate  the  use  of  the  vectors  and  the  algorithm  used  to  develop  the 
signal 

=  «0  Kt  -  To)  +  CX^s(t  -T^)  +  a2S(t-T2)  +  a^t  -  Tj) 

where 

r  =  [l  10  12  20] 
a  =  [l  0.5  0.3  O.l] 

thus 

=  s(t)  +  0.5- s(t  -9-dt)  +  03^s{t  - 1 1  ^dt)  +  0.1  •  5(?  - 19  •  dt) 
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Figure  10  Multipath  Example  Signal 

The  top  graph  in  Figure  10  shows  the  undelayed  signal,  the  middle  graph  shows 
the  direct  path  and  multipath  signals  plotted  individually  and  the  bottom  graph  shows  the 
sum  of  the  individual  and  multipath  signals. 

If  the  multipath  delay  time  and  amplitude  are  increased  with  the  following 
parameters 
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combined  composite  undelayed 


r  =  [l  20  22  30] 
a  =  [l  0.5  0.3  0.5] 

the  resulting  multipath  signal  will  be  as  shown  here: 


Multipath  simulation 


Figure  1 1  Multipath  Example  with  Large  Delay 

Notice  that  in  both  examples  the  multipath  signal  has  a  different  amplitude  from 
the  original  signal  and  that  the  larger  the  delay  the  larger  the  change  in  the  phase  of  the 
signal.  Below  is  an  example  of  destructive  interference  obtained  with  the  values 
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Figure  12  Multipath  Example  with  Destructive  Interference 


This  chapter  has  described  the  multipath  environment  and  presented  a  simulation 
model  to  generate  a  signal  based  on  user  defined  parameters.  The  next  chapter  will 
describe  a  Kalman  filter  which  is  used  for  tracking  the  GPS  signal. 
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IV.  KALMAN  FILTER  FOR  TRACKING  PHASE 


A.  DISCRETE  KALMAN  ESTIMATOR 

The  development  of  a  filter  to  track  the  GPS  signal  in  a  multipath  environment 
began  with  the  implementation  of  a  Discrete  Kalman  Estimator  for  tracking  the  phase 
changes  of  BPSK  signal.  The  derivation  below  closely  follows  that  given  by  Professor 
Kirk  (1975). 


1.  Plant  Equations 


The  plant  equations  were  derived  as  follows: 
X,  =  cos(fl?  t) 

=x,  =-G)sin(6)0 

f X,'\  (  COS(Q}t)  ^ 

x(0=  ’  =  \ 


-msin(oO^_r  ^  f  cos(<nO  ^ 
cos{cot))  oj  V-<ysin(ftjO/ 


2.  Plant  Characterization 


The  discrete  model  of  the  plant  is  characterized  by  the  linear  discrete  difference 


equations: 


x(k  + 1)  =  ^  :r(A:) + h'(A:) 


2(k)  =  Cx(k)+v(k) 

where: 

x(k)  is  the  n-dimensional  state  vector  at  time  k, 

z(k)  is  a  q-dimensional  output  measurement  vector  at  time  k 

w(k)  is  a  p-dimensional  vector  of  random  forcing  inputs  at  time  k 


(t)  is  the  state  transition  matrix 
C  is  the  observation  matrix 
and 

v(k)  is  the  q-dimensional  vector  of  random  measurement  noise  at  time  k 

3.  Gain,  and  Covariance  Equations 

The  gain,  G,  and  covariance  equations  derived  by  Kirk  (1975)  are: 

Gik)  =  P{k\k  - 1)  •  •  [CPik\k- 1)  •  +  rY 

P(k\k)  =  [/  -  Gik)  •  C] •  Pik\k  - 1) 

P(k+l\k)  =  ^-Pik\k)-^^+0 

where 

C  is  the  observation  matrix 
I  is  the  identity  matrix, 

G  is  the  Kalman  gain 
P  is  the  variance  of  estimation  error 
Q  is  random  process  forcing  input  covariance  matrix 
R  is  the  variance  of  measurement  noise. 

The  estimator  equations  are:  xik\k)  =  xik\k  - 1)  +  Gik) •  [zik)  -  C •  xik\k  - 1)] 
xik\k-\)=^<f>-xik-\\k-\)  +  K‘Uik-\) 
with  the  initial  condition  x(0|-l)  =  x^ 


^  ^ij\j )  ^  j  given  j  and  is  defined  as  the  estimate  of  x  at  time  j  given  measurements  up  to  and 

including  time  j. 
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4.  Assumptions 


a)  The  measurement  noise  has  zero  mean; 

£[v(yfc)]  =  0  A:  =  0,1,...^ 

is  uncorrelated  and  has  covariance  R(k): 

£[v(A:)v^(7)]  =  £[v(7)v^(A:)]  =  R{k)5^ 

where  5 is  the  Kronecker  delta  fimction  defined  by 

^  ^  =  ^ 

b)  The  initial  state  is  a  random  variable  with  known  mean 

E[x(o)]=;; 

and  covariance; 

£j^[x(0)-^][x(0)-^;f]  =  M 

The  measurement  noise  and  initial  state  are  uncorrelated: 

£[x(0)v'’(yt)]  =  £[v(A:)x'’(0)]  =  0  k  =  OX 

c)  by  the  linear  relationship 

x(A^|  A:)  =  x(Arl  A:  - 1) + G(A:)  •  [z(A^)  -  C  ■  x{k\k  - 1)] 

where  G(k)  is  the  Kalman  gain  and  z(k)  -  C  ■  x(A:|A:  - 1)  is  the  correction  term,  or 

the  difference  between  the  measured  and  predicted  value. 

d)  The  random  process  forcing  input  has  zero  mean 

£'[w(A:)]  =0  ^  =  0,1, . . . 

is  uncorrelated  and  has  covariance  Q(k) 

£[^h-(A:)w^(7)]  =  £[w(7)w^(A:)]  =  Qik)S^  j,k  =  0,1,2,... 

e)  The  forcing  random  process  and  initial  state  value  are  uncorrelated 

^  Ell  is  Ihe  expected  value  or  mean  of  the  random  variable. 
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E\wik)x^  (0)]  =  £:[x(0)w'’(^)]  =  0  k  =  0,1,2, . . , 
f)  The  forcing  random  process  and  measurement  noise  process  are  uncorrelated 
E\w{ky  (y)]  =  E\v(J)w^  (A:)]  =  0  k  =  0,1,2,  . 


B.  COMPUTATIONAL  PROCEDURE 

The  computational  procedure  is  shown  in  the  following  algorithm. 

Initialize  with 
k=0 

x(ohi)=^: 

P{k\k-\)  =  M 

The  procedure  iterates  on  k  in  the  loop: 

G{k)  =  P{k\k-\)  C^  [C- P{k\ k-\)-C^+  i?]"’ 

P{k\k)  =  [/ - G(k)  ClPmk-l) 

P{k  + 1|^)  =  P(k\k)  -^^+0  {this  will  be  P{k\k  - 1)  next  time  through) 

take  measurement  z(t) 

x{k\k)  =  x{k\k  - 1)  +  G{k)  •  [z{k)  -  C  •  x(k\k  - 1)] 
x(k\k  -\)=<j>- x(k  -  \\k  - 1)  +  ^-u(k  - 1) 
k=k+l 

End  of  the  Loop 


C.  PHASE  TERMS 


The  QPSK  signal  can  be  represented  by 


28 


s(t)  =  A  cos(fi>r  +  ^  (0) 

where  ^(t)  represents  the  phase  change  of  the  signal 


<l>o 

2 

<l)^+n 

A 

+ — 


Substituting  in  these  phase  values  and  simplifying  using  trigonometric  identities 
results  in 


Acos(o}t  +  ^o) 

-Asin(cot  +  ^J 

5(0  =  1  .  ,  ,  . 

-Acos(cot  +  <l>^) 

Asva{o)t-\-^o) 

Using  these  relationships  and  the  sine  and  cosine  information  in  the  state  estimates, 
the  Kalman  filter  can  estimate  which  phase  change  has  occurred.  The  error  term  is 
calculated  by  finding  the  estimate  of  the  signal  that  minimizes  the  difference  between  the 


measured  value.  Assuming 


i(A:|^-l)  = 

C  =  [1,0] 
C,.=[0,1] 


cos(®/  +  ^^)  ^ 
-wsva.i(ot  +  ^g)) 


the  four  error  terms  can  be  calculated  as  shown  below,  and  the  minimum  value  is 
used  to  update  the  prediction 

=  z{k)  -  Cx(A:|  - 1)  =  z(k)  -  cos(«u  A  +  ) 

=  z{k)  +  Cp,x(A|  A  - 1)  /  fi)  =  z{k)  -  sin((y  A  +  (p^) 

=  z(k)  +  Cx(A|  A  - 1)  =  2(A)  +  cos(fi)  A  + 

£3  =  2(A)  -  C3„x(A1  A  - 1)  /  ®  =  2(A)  +  sin((a  k  +  (f>^) 
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These  error  terms  are  then  used  to  update  the  estimate  of  the  signal  without  any 
phase  shifting  which  will  be  called  y,  and  an  estimate  of  the  phase  shifted  signal  called 
y_ps. 

Next,  the  results  of  using  the  developed  Kalman  filter  to  track  the  phase  of  a 
BPSK,  QPSK,  and  unbalanced  QPSK  signal  will  be  shown. 
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V.  INITIAL  RESULTS 


This  chapter  will  show  the  results  obtained  when  the  Kalman  filter,  discussed  in  the 
previous  chapter,  is  used  the  track  the  phase  of  a  BPSK,  QPSK  and  unbalanced  QPSK^ 
signal.  The  performance  of  the  Kalman  filter,  will  be  shown  with  and  without  noise,  and 
in  the  presence  of  multipath  interference.  It  will  be  shown  that  the  Kalman  filter  will  track 
the  combined  line  of  sight  and  direct  path  signals  phase  and  amplitude.  It  will  then  be 
shown  that  when  given  a  priori  information  the  filter  can  distinguish  between  the  line  of 
sight  and  multipath  signals  and  correctly  track  them,  given  a  static  multipath  environment. 


A.  SIGNAL  GENERATION 

The  Matlab  m-file  uij>ars.m  sets  up  the  user  interface  shown  in  Figure  13  that  can 
be  used  to  control  the  simulation. 

Each  entry  can  be  edited  by  the  user  to  change  each  of  the  parameters  discussed 

below: 

•  Signal  freq  -  Sets  the  frequency,  of  the  oscillator,  in  Hertz,  used  by  the 
carriers  in  the  BPSK  and  QPSK  signal  generation. 

•  Phase  Offset  -  Determines  the  initial  phase  offset,  of  the  line  of  sight  signal  in 
degrees. 

•  Amplitude  -  Sets  the  amplitude  of  the  generated  signal. 

•  -  This  sets  the  sampling  period  in  seconds.  In  the  example  below  Ts  =0.01, 
resulting  in  100  samples  per  second. 

•  Stop  time  -  Sets  the  Simulink  parameter  stop  time  to  control  length  of 
simulation. 

•  Tb  -  Sets  the  bit  period  which  determines  the  length  of  time  for  zero-order  hold 
per  data  bit. 


''  An  unbalanced  QPSK  signal  is  a  QPSK  signal  with  differing  energy  in  the  in-phase  and  quadrature 
phase  components.  For  GPS  the  LI  signal  the  amplitude  of  the  quadrature  term  is  equal  to  times  the 
in  phase  term. 
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•  SNR  dB  -  Using  this  value  and  the  entered  amplitude  value  of  the  signal,  the 
random  number  generator  is  multiplied  by  a  scaling  factor  to  simulate  the 
desired  SNR  level. 

•  tau  -  Delay  vector  used  to  control  the  multipath  delay  times.  See  Chapter  HI 
Section  C  for  more  information. 

•  alpha  -  Amplitude  vector  -  Used  to  control  the  multipath  amplitudes. 


Figure  13  User  Interface  to  Control  BPSK  and  QPSK  Simulation 
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The  pulldown  menu  Init  initializes  all  the  entered  values.  Plot  gives  the  user  a 
choice  of  plots  to  choose  from  to  document  the  output  of  the  simulation.  Sim  allows  the 
user  to  choose  which  simulation  will  be  run. 

B.  BPSK 

Using  Matlab  and  Simulink,  a  BPSK  signal  was  generated  and  the  algorithm 
discussed  in  the  previous  chapter  was  used  to  track  the  phase  changes  in  the  signal.® 

1.  BPSK  Signal  Without  Noise 

Figure  14  shows  the  performance  of  the  algorithm  when  there  is  no  noise  present. 
In  this  graph,  x  is  the  plant  state,  z  is  the  measured  signal  (x  +noise),  y  is  the  recovered 
carrier,  and  y_ps  is  the  recovered  phase  shifted  signal.  The  last  graph  shows  the  error  in 
the  recovered  signal  which  is  the  difference  in  the  generated  signal  and  recovered  estimate, 
x-y_ps.  Notice  there  is  an  initial  error  during  initialization,  however  the  error  quickly  goes 
to  zero  as  expected  in  the  absence  of  noise. 

2.  With  Noise 

Keeping  all  variables  the  same  as  above  and  adding  zero  mean  Gaussian  white 
noise,  by  using  the  random  function  in  Matlab,  with  an  SNR  value  of  12  dB®,  the  carrier  is 
still  recovered,  however  small  errors  are  now  introduced  in  the  estimate  of  the  phase 
shifted  signal. 


^  See  Matlab  m-files  kaljnit.m,  kalman2.m  and  the  Simulink  program  kal  sim 
®  Spilker  states  that  typical  values  for  C/No  are  40.6  dB-Hz,  which  corresponds  to  an  SNR  of  25  dB.  See 
Appendix  C  for  a  link  budget  estimate  of  the  signal  to  noise  ratio. 
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Kalman  Filter 
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Figure  14  Kalman  Filter  as  Phase  Lock  Loop  for  BPSK  Without  Noise 
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time  -  sec 

Figure  15  Kalman  Filter  for  BPSK  with  Noise 


C.  QPSK 

Next  the  results  for  a  QPSK  simulation  are  shown.  These  plots  were  generated 
using  kalqsim.m 

1.  Without  Noise 

Figure  16  shows  the  performance  of  the  algorithm  when  there  is  no  noise  present. 
The  true  signal,  x,  and  the  estimate  y_ps  and  the  error  in  the  estimate,  x-y_ps  are  all 
plotted. 
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Kalman  Filter  Error 


Figure  16  Kalman  Filter  for  QPSK,  Without  Noise 


2.  With  Noise 

Figure  17  shows  the  estimation  when  white  noise’  is  added  with  an  SNR  value  of 
12  dB.  Figure  18  shows  the  error  in  the  estimation.  The  carrier  is  still  recovered,  however 
small  errors  are  now  introduced  in  the  estimate  of  the  phase  shifted  signal.  The  cause  of 
the  error  will  be  explained  in  the  next  section. 


’  Throughout  this  stuffy  only  gaussian  white  noise  was  considered.  Colored  noise  from  other  sources,  such 
as  electromagnetic  interference  was  not  considered.  It  was  assumed  that  any  such  noise  would  have  the 
same  affect  at  both  receivers  and  would  not  change  the  difference  in  the  measured  phase  values. 
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D.  SOURCES  OF  ERROR  IN  NOISY  CASE 


The  error  in  the  estimation  caused  by  the  noise  can  best  be  explained  by  examining 
the  BPSK  simulation.  In  the  algorithm,  the  Kalman  filter  tracks  the  carrier  signal  and 
estimates  which  generated  phase  shifted  version  has  the  minimum  error  when  compared 
with  the  measured  signal  z.  Figure  19  is  a  plot  of  the  original  BPSK  signal,  x,  the  signal 
with  noise  added,  z;  the  estimated  signal,  y_ps;  the  two  possible  signals,  xO  and  x2,  that 
must  be  decided  between;  and  the  error,  eO  and  e2,  used  as  the  basis  for  the  decision. 

The  error  in  the  signal  is  best  described  as  “clicks”  which  occur  when  the  Kalman 
filter  is  unable  to  distinguish  between  the  correct  and  incorrect  phase  because  they  both 
have  approximately  the  same  value.  The  reader  is  invited  to  examine  this  effect  between 
times  1  and  2  seconds.  The  correct  estimated  signal  is  X2  which  results  in  the  minimum 
error  over  the  bit  period,  however  at  approximately  time  1 .9  the  out  of  phase  signal 
estimate  XO  has  the  same  value  as  X2  and  the  error  eO  is  less  than  el  due  to  the  added 
noise.  This  produces  an  erroneous  decision  to  switch  phase.  Typically  the  erroneous 
decisions  all  occur  when  both  the  correct  and  incorrect  signal  estimates  are  near  a  zero 
crossing.  This  is  due  to  the  nature  of  the  cosine  fimction  which  results  in 
cos(90°)=cos(270°)=0. 
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E.  UNBALANCED  QPSK 


As  discussed  in  Chapter  II,  the  GPS  LI  signal  is  actually  an  unbalanced  QPSK 
signal  with  dififering  energy  in  the  quadrature  and  in-phase  components.  As  shown  in 
Chapter  11  Section  F  on  page  15  the  LI  signal  is  given  by 

Lj  (ft)j  t)  =  [P(/)  0  D(0]  cos(fiJi  t)  -I-  [G(/)  0  iX^)]  sin(fi?,  t) 

An  example  of  the  unbalanced  QPSK  signal  is  shown  in  Figure  20.  The  top  graph 
shows  the  composite  of  the  binary  data  representing  [G(/)  0Z)(Q]  with  the  quadrature 

component,  the  middle  plot  shows  the  binary  data  representing  [P(Q0Z)(Q]  with  the  in- 

phase  component.  In  both  cases  the  Ai  has  been  normalized  to  1 .  The  last  plot  shows  the 
sum  of  the  two  signals  which  is  the  resulting  unbalanced  QPSK  signal. 


QPSK  Simulation 


Figure  20  Unbalanced  QPSK  Signal  Generation 
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When  the  Kalman  filter  is  used  to  track  this  signal  as  shown  in  Figure  21,  errors 
occur  when  the  in-phase  and  quadrature  terms  have  differing  data.  The  filter  locks  on  to 
the  initial  ampUtude  and  assumes  the  amplitude  will  not  change.  In  this  example  the 
absolute  value  of  the  amplitude  is  the  same  when  the  data  bits  are  the  same  (both  +1  or 
both  -1)  and  is  different  with  the  data  bit  have  opposite  values.  The  state  estimation  was 
developed  under  the  assumption  the  signal  would  have  a  constant  amplitude.  In  order  to 
resolve  this  problem  the  logic  for  phase  change  detection  would  need  to  be  modified. 


Kalman  Filter 


Figure  21  UQPSK  Without  Noise 
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When  noise  is  added  to  the  unbalanced  QPSK  as  shown  in  Figure  22  the 
perfonnance  is  the  same  as  in  the  QPSK  example  when  the  data  bit  for  the  in-phase  and 
quadrature  phase  term  is  the  same,  however  when  they  are  different,  the  performance 
degrades.  The  error  in  the  estimation  is  shown  in  Figure  23. 


Kalman  Filter 
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Figure  22  UQPSK  with  Noise 
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Kalman  Filter  Error 


Figure  23  Error  in  UQPSK  Estimator 


F.  MULTIPATH  RESULTS 


1.  Without  A  Priori  Knowledge 

The  multipath  parameters  used  to  generate  Figure  10  are  used  to  generate  the 
measured  signal  z  and  the  results  obtained  are  shown  below  in  Figure  24.  The  estimate, 
yjps,  and  the  direct  path  signal,  x,  are  shown  plotted  together  and  the  error,  or  the 
difference  between  x  and  y_ps  is  also  shown. 
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Kalman  Filter  Error 


time  ~  sec 


Figure  24  Results  of  Kalman  Filter  Estimation  in  Presence  of  Multipath 


The  Kalman  filter  has  correctly  determined  the  phase  changes,  however  due  to  the 
multipath,  the  amplitude  of  the  estimation  is  now  different  from  the  original  signal.  Using 
the  larger  multipath  spread  parameters  used  to  generate  Figure  1 1  the  resulting  phase 
error  is  observed  as  shown  Figure  25. 
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Kalman  Filter  Error 


Figure  25  Results  of  Kalman  Filter  Estimation  in  Presence  of  Multipath  with  Increased 

Delay 


The  results  of  these  last  two  plots  indicate  the  Kalman  filter  performance  in  the 
presence  of  multipath  is  much  worse  than  in  the  presence  of  white  noise.  If  multipath  is 
present  it  will  introduce  a  phase  error  that  will  be  unique  to  each  antenna.  When  the 
differential  phase  is  calculated  it  will  now  include  the  multipath  phase  error  and  will  no 
longer  be  just  a  function  of  the  geometry. 
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2. 


With  A  Priori  Knowledge 


If  the  multipath  amplitude  and  time  delay  are  known,  the  Kalman  filter  can  be 
modified  to  track  both  the  the  direct  path  signal  and  delayed  multipath  signals.  A  simple 
case  with  a  single  multipath  channel  is  shown  below  with  the  following  assumptions: 

The  combined  signal  is  given  by: 

The  incoming  signal  is  not  changing  phase  and  is  given  by 

r(#)  =  A  cos(o  t  +  $it)) 

where  6(t)=  a(0f  and  a(t)  =  0,1, 2, or  3  representing  the  intelligence  for  the  BPSK 
or  QPSK  modulation. 

The  time  delay  r  is  converted  into  the  phase  delay 

s{t  -t)  =  A  cos((v(f  -t)  +  6{t  -t))  =  A  cos(cd  i-coT  +  0(^-r))  =  A  cos(co  t  + 
and  the  composite  multipath  signal  is  now  given  by 

s^(t)  =  Acos{Q)t  +  0)-\-  a  Acos{(ot  +  <!>) 

The  plant  equations  are  now 

-  A  cos(o  t  +  ff) 

X2  =  =  -CO  A  sin(ta  t  +  0) 

X3  =  a  A  cos{co  t  +  ^) 

X4  =  Xj  =-coaA  sin(®  t  +  ^) 


x(0  = 


hi 

^  A  cos(o)  t  +  0)  ^ 

^2 

-CO  A  sin(o  t  +  0) 

— 

a  A  cos(co  t  +  ^) 

U4J 

^-coaA  sin(fi>  t  +  ^)J 

^  -6)  Asin(a)t  +  9)  '' 

0 

0 

^  A  cos(e?  t  +  0)  ^ 

X2 

-  co^A  cos(a)  t0) 

-co^  0  0  0 

-  (oA  sin(o  t  +  0) 

-coa  A%m{0t +  (!>) 

0  0  0  1 
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U4J 

K-co^aA  cos(<o  t  +  ^)) 

V  0  0  -  Oy 

^-coa  Asin(a)t +  ^)j 

=  Ax(t) 
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Since  a  priori  knowledge  is  known  the  values  for  the  estimate  at  time  zero  can 

now  be  initialized  as  follows: 

A 
0 

a  A  cos(^) 

-oa  A^n{<f>) 

Figure  26  shows  the  line  of  sight  signal,  delayed  signal,  and  the  measured  signal  z 
which  were  generated  using  a=0. 5  and  t=10  with  Ts=.  1 .  For  this  example  the  multipath 
environment  is  assumed  to  be  static.  In  reality,  the  multipath  environment  would  be  much 
more  dynamic  due  to  the  satellites  orbit,  the  satellites  changing  attitude,  and  the  changing 
position  of  the  GPS  satellites. 


Figure  26  Example  Signal  for  A  Priori  Case 


A1 


los  &  xkl<ps(1 .:)  (estimate  = 


Figure  27  A  Priori  Knowledge  Results  Showing  Signals  and  Estimates 


Figure  27  shows  the  estimations  and  original  signals  plotted  on  top  of  each  other. 
This  demonstrates  how,  given  a  priori  knowledge,  the  line  of  sight  and  delayed  signal  can 
be  tracked.  Now  there  is  no  phase  error  due  to  the  multipath  signal,  unlike  the  examples 
shown  without  a  priori  knowledge  when  the  filter  tracked  the  phase  of  the  combined 
signal. 

When  noise  is  added  in  this  simple  case  with  12  dB  SNR  as  shown  in  Figure  28 
and  Figure  29  there  is  also  improved  performance  with  a  priori  knowledge. 
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Figure  28  Generated  Signal  for  A  Priori  Example  with  Noise 
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Figure  29  A  Priori  Knowledge  with  noise  Results 


G.  EFFECT  OF  Q  MATRIX  ON  KALMAN  GAINS 


The  random  process  forcing  input  covariance  matrix  Q  is  a  measure  of  the 
expected  excitation  of  the  noise.  It  is  used  in  the  calculation  of  the  Kalman  gains  as 
shown  in  Chapter  IV,  Section  A-3  on  page  26  which  are  restated  here; 

G{k)  =  Pik\k  - 1)  •  C""  •  [C  •  P{k\k  - 1)  •  C  +  r]' 

P(k\k)  =  [I-G(k)C]P(k\k-l) 

P(k  +  l\k)  =  (^-P(k\ky,^^+Q 


The  Q  matrix  is  a  function  of  the  state  noise  coefiScient  matrix,  F ,  and  the 
variance  of  the  forcing  function,  w,  as  shown  below. 

dt^ 


r  = 


Q 


Q  = 


2 

dt 

r-r^ 

dt'^ 

IT 

dt^ 


var[>v] 

dt^ 


2 

dt^ 


var[w] 


As  there  is  no  forcing  function,  w,  the  value  of  Q  was  initially  set  to  zero,  however 
in  order  to  determine  its  impact,  different  values  for  the  variance  of  w,  var[w],  were  used. 
The  filter  results  and  Kalman  gain  values  G1  and  G2  are  shown  in  Figure  30  through 
Figure  33  for  different  values  of  var[w]  and  Q.  As  the  value  for  var[w]  is  increased,  the 
steady  state  gain  increases.  The  only  noticeable  effect  occurs  when  var[w]  is  greater  than 
100.  With  an  increase  in  var[w]  and  Q,  there  is  also  an  increase  in  the  Kalman  gains.  This 
results  in  more  responsiveness.  If  the  gain  values  are  too  large  the  filter  will  be  too 
responsive  resulting  in  more  errors. 
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Kalman  Filter  -  Signal  Options 


Figure  30  Effect  of  Gain  Values  with  var[wj=0,  Q 
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Figure  34  Steady  State  Kalman  Gains  with  var[w]=150 


The  steady  state  gain  values  after  time  t=2  seconds  are  shown  in  Figure  34.  For  all 
other  values  of  var[w]  shown,  the  corresponding  steady  state  Gain  values  approached 
zero.  ’ 


H.  SUMMARY 

This  chapter  has  presented  the  results  of  using  the  Kalman  filter  to  estimate  the 
phase  in  BPSK,  QPSK,  and  Unbalanced  QPSK  signals  to  simulate  the  various  formats  of  a 
GPS  signal.  Results  have  been  shown  for  the  performance  with  and  without  Gaussian 
white  noise  and  with  multipath  interference. 

The  designed  filter  works  best  with  the  BPSK  signal.  Gaussian  white  noise 
degrades  performance  due  to  the  effects  of  the  correct  and  incorrect  phase  having  the 
same  value  at  a  zero  crossing.  It  may  be  possible  to  modify  the  phase  switching  logic  to 
reduce  this  effect  by  allowing  a  phase  change  to  only  occur  at  or  near  the  expected  time  of 
the  end  of  a  bit  period  or  when  the  error  differences  are  greater  than  some  minimum  value. 
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The  multipath  interference  was  shown  to  result  in  errors  in  the  estimate’s  phase 
and  amplitude.  This  error  is  caused  by  the  filter  tracking  the  phase  of  the  combined 
multipath  signal.  It  was  shown  that  with  a  priori  knowledge  of  multipath  delay  and 
amplitude,  the  filter  is  able  to  discriminate  between  the  direct  path  and  the  delayed  signals 
and  to  eliminate  this  error.  Sources  of  the  a  priori  information  might  include  ground 
testing  by  placing  a  receiver  into  a  self-survey  mode  and  allowing  the  GPS  satellites  to 
pass  overhead.  (Lightsey)  It  may  also  be  possible  to  use  the  technique  suggested  by 
Kumar  and  Lau  to  estimate  the  multipath  parameters  through  seismic  deconvolution 

For  this  case  it  was  assumed  that  the  multipath  environment  was  not  changing.  In 
reality  the  multipath  environment  would  be  very  dynamic  as  a  result  of  the  satellites  orbit 
and  changing  attitude.  Additionally,  the  attitude  determination  and  control  system  only 
uses  4  out  of  possibly  1 1  visible  satellites.  The  satellites  selected  would  also  effect  the 
multipath  environment.  To  overcome  this  may  require  frequent  re-initialization  of  the 
Kalman  filter. 
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VI.  CONCLUSION 


A.  SUMMARY 

This  study  has  focused  on  the  carrier  phase  signal  which  is  used  by  GPS  based 
attitude  detemaination  systems.  A  GPS  carrier  signal  and  a  multipath  environment  were 
described  and  modeled,  A  Kalman  filter  was  developed  to  track  the  GPS  signal  and 
results  were  shown  for  Gaussian  noise  and  in  the  presence  of  multipath. 

The  developed  Kalman  filter  was  shown  to  work  best  for  the  L2,  BPSK  signal  due 
to  the  minimum  complexity  involved  in  the  phase  changing  logic.  It  was  shown  that  the 
filter  is  capable  of  discriminating  between  the  line  of  sight  and  delayed  signals  and 
performed  well  in  the  presence  of  multipath,  given  good  a  priori  knowledge. 

More  refinement  is  required  in  the  phase  change  detection  logic  to  eliminate  errors 
caused  by  the  error  terms  having  the  same  value  at  a  zero  crossing.  The  Q  excitation 
matrix  needs  considerable  adjustment  to  handle  changing  phase  in  the  presence  of 
multipath.  The  source  of  the  a  priori  information  needs  to  be  developed. 

B.  OTHER  METHODS  OF  MULTIPATH  MITIGATION  FOR  POSSIBLE 

STUDY 

As  multipath  accounts  for  90%  of  the  error  in  an  attitude  determination  system, 
much  work  is  currently  being  done  in  this  area.  Three  categories  of  multipath  mitigation 
methods  that  have  been  used  include:  1)  Methods  involving  the  geometry  of  the 
multipath,  2)  Analyzing  the  consistence  of  the  different  code  measurements  and  3) 
Exploiting  information  coming  fi’om  different  chaimels. 

Of  these  three,  the  easiest  to  implement  involves  geometry.  One  approach  is  to 
focus  on  improving  the  antenna  or  its  position.  This  may  be  accomplished  by  placing  the 
antenna  in  isolated  locations  minimizing  the  multipath,  something  that  is  not  always 
possible  on  a  satellite.  Another  method  of  multipath  rejection  is  achieved  by  reducing  the 
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left-hand  circularly  polarized  gain  (LHCP)  of  the  antenna  without  reducing  the  right-hand 
circularly  polarized  gain  (RHCP).  This  is  because  the  satellite  signals  are  RHCP  and  the 
reflected  multipath  signals  tend  to  be  either  linearly  polarized  or  even  LHCP,  depending 
on  the  reflecting  surface.  (Van  Dierendonck)  An  example  of  an  antenna  designed 
specifically  to  reduce  multipath  effects  is  the  choke  ring  antenna.  This  anteima  is 
comprised  of  vertical-aligned  concentric  rings  centered  about  the  antenna  element  that  are 
cormected  to  the  ground  plane.  The  vertical  rings  shape  the  antenna  pattern  such  that  the 
multipath  signals  incident  on  the  antenna  at  the  horizon  and  negative  elevation  angles  are 
attenuated  (Kaplan). 

Other  methods  of  dealing  with  the  multipath  problem  involve  modification  to  the 
GPS  receiver  for  code  measurements.  One  approach  involves  reducing  the  early-late 
delay  spacing  among  the  correlators  in  the  GPS  receiver  code  lock  loop.  This  technique 
provides  limited  results  if  the  early-late  spacing  is  smaller  than  the  initial  delay  error  due  to 
multipath.  (Kumar) 

An  example  of  the  third  category  is  the  method  presented  by  Kumar  and  Lau 
involving  the  use  of  optimal  deconvolution  to  estimate  the  impulse  response  of  the 
effective  multipath  chaimel  and  obtain  an  inverse  filter  which  equalizes  the  multipath 
channel  (Kumar)  Another  technique  involves  adaptively  estimating  the  spectral 
parameters  (frequency,  amplitude,  phase  offset)  of  multipath  in  the  associated  SNR  and 
then  constructing  a  profile  of  the  multipath  error  in  the  carrier  phase.  A  multipath 
correction  is  then  made  by  subtracting  the  profile  from  the  actual  phase  measurement. 
(Comp) 
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APPENDIX  A  MATLAB  PROGRAMS 


A.  DO_KAL.M 


%  File  Name:  do_kal.m 
%  Last  Updated:  18  Apr  97 
%  Written  By:  T.  H.  Newman 
%  Operating  Environment:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description:  Script  file  to  set  up  and  run  kalman  filter  simulations 

clear 

ui_vars 

%  End  of  File  do_kal.m 


B.  G1.M 


%  File  Name:  gl.m 
%  Last  Updated:  18  Oct  96 
%  Written  B>’:  T.  H.  Newman 
%  Operating  Environment:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description:  This  function  generates  the  Gl(t)  signal 
%  for  GPS  C/A  Code  generator 
% 

% 


function  [new_R,  g]  =  gl(R,  xl) 
if  xl==l 

R  =  ones(size(R));  %  Reset  to  all  ones  at  XI  epoch 
end; 

new_R=gen_poly([l,3,10],R); 

g=R(10); 


%  End  of  File  gl.m 
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c. 


G2.M 


%  File  Name:  g2.m 
%  Last  Updated:  18  Oct  96 
%  Written  By:  T.  H.  Newman 
%  Operating  Enviroiunent:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description;  This  function  generates  the  G2(t)  signal 
%  for  GPS  C/A  Code  generator 
% 

% 


function  [new_R,  g]  =  g2(R,  xl,  sv) 


ifxl==l 

R  =  ones(size(R));  %  Reset  to  all  ones  at  XI  epoch 
end; 

new  R  =  gen_poly([l  2  3  6  8  9  10],  R); 

%new_R(l)  =  mod2(R(2)+R(3)+R(6)+R(8)+R(9)+R(10)); 
%new_R(2;10)  =  R(l:9); 


%  phase  selert  logic,  see  p.  92  KAPLAN 
%  This  changes  based  on  satellite  choosen. 
sv 'tap  =[26 
37 
48 
59 
19 
2  10 
1  8 
29 
3  10 
23 
34 
56 
67 
78 

8  9;  9  10;1  4;2  5;3  6;  4  7;  5  8;  6  9;  1  3;  4  6;  5  7;  6  8;  7  9; 
8  10;  1  6;  2  7;  3  8;  4  9;  5  10;  4  10;  1  7;  2  8;  4  10]; 


g=xor(  R(sv'_tap(sv,  1)  ),  R(sv_tap(sv,2))  ); 


%  End  of  File  g2.m 
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D. 


GEN  POLY.M 


%  File  Name;  gen_poly.m 
%  Last  Updated:  18  Oct  96 
%  Written  B)':  T.  H.  Newman 
%  Curating  Environment;  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description:  This  function  generates  the  linear  code  generators 
%  the  polynomial  is  of  the  form  1  +  sum(x(i))  where 
%  sum(x(i))  means  that  the  ou^t  of  the  i'th  cell  of  the  register 
%  is  used  as  the  input  to  the  m^ulo-2  adder  (exclusive  or) 

%  and  the  1  means  that  the  output  of  the  adder  is  fed  to  the  first  cell 
% 

%  so  1  +  xl  +  x3  +  xlO  could  be  generated  with  poly  =[1  1  3  10]  and  R 
%  is  the  register  to  perform  the  operation  on 

function  r  =  gen_poly(poly,  R) 

result  =  mod2(  sum(R(poly(2:length(poly) ) )));  %  xor  registers  according  to 

%  polynomial 

r(2:length(R))  =  R(l:Oength(R)-l));  %  shift  registers  right 

r(poly(l))=result;  %  put  result  in  register  according 

%  to  polynomial 


%  End  of  File  gen_poly.m 


E.  INTT^I 


%  File  Name:  init.m 
%  Last  Updated;  11  Nov  96 
%  Written  By:  T.  H.  Newman 
%  Operating  Enviromnent:  Sun  Qpenwin  &  Win  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description:  This  script  file  sets  the  parameters  for 
%  the  GPS  signal  simulation 

global  w; 

%  my  _glb; 
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fo  =  10.22999999543e6; 

fl  =  120*fo;  %  frequency  of  L2  signal 

f2  =  154*fo;  %  frequency  of  LI  signal 


To  =  .001; 

codejength  =  10; 

[s  1  ,s2.s31=sig  genfcode  length): 
sigl=[(l;length(sl))',sr]; 
sig2=[(l:length(s2))',s2*]; 
sig3=[(l;length(s3))',s3’]; 

kal_init; 

%  End  of  File  init.m 


F.  KAL  INIT.M 


%  File  Name;  kal_init.m 
%  Last  Updated:  27  Feb  97 
%  Written  By:  T.  H.  Newman 
%  Operating  Environment:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description;  Initialization  file  for  Kalman  Filter 
% 


global  Pkkml; 
global  R; 
global  Q; 
global  I; 
global  Phi; 
global  xkkml; 
global  Tb; 
global  Noise_gain; 
global  dt; 

global  w; 

%  myglb; 

A=[0  1;  -vv^2  0];  %  xl  =  cos(wt) 

%  x2  =  xldot  =  -w*sin(wt) 
%  x2_dot  =  -w'^2*cos(wt) 

%  X  =  [xl  x2]' 
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B=[0  w'^2]'; 


%  x_dot  =  A*[xl_dot,  x2_dot] 
%  =>A=[0  l;-w''2  0]’ 

%  x_dot  =  Ax  +  Bu 
%  u  is  control  input 


C=[l  0];  %  will  use  to  get  xl  estimate 

D=0; 


I=ej'e(2);  %  Identity  Matrix 

[Phi,Del]=c2d(A,B,dt);  %  converts  the  continuous-time  system; 

%  x_dot  =  Ax  +  Bu 
%  to  the  discrete-time  state-space 
%  system: 

%  x[n+l]  =  Phi  *  x[n]  +  Del  *  u[n] 
%  assuming  a  zero-order  hold  on  the 
%  inputs  and  sample  time  dt 


Pkkml=le6*I; 


R=l; 


%Q=0*I; 


%  Initialize  P(k|k-1)  =M 

%  where  M  =  covariance  of  initial  state 
%  want  this  number  large  to  give 
%  weight  to  1st  measurment 
%  var[e(010)]  =  var[v(0)]  =  R 
%  V  is  reference  input 
%  (noise  in  measurement) 

%  Q  is  measure  of  exitation  by  forcing  function 


var_w=150; 

Q  =  [  (dt'^4)/4  (dt'^3)/2;  (dt^3)/2  dt^2  ]*var_w; 


xkk  =  [0  0;  0  0] ;  %  initialize  x(k|k)  to  all  zeros 


xkkps  =  [0  0;  0  0]; 


% 


x(k|k)  is  state  of  plant  at  time  k  given  k  observ'ations 


xkkml  =  [0,0;  0,0]; 


%  initialize  x(k|k-l)  =>  x(0|-l)  =  x_o  =  0 


%  clear  all  the  variables  used  for  saving  plotting  information 

clear  my_rand; 

clear x; 

clear  sm; 

clear  xdelay; 

clear  z; 

clear  y; 

clear  xkk_ps; 

clear  t; 
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%EndofFilekal  initm 


G.  KALMAN.M 


%  File  Name:  kalman.m 
%  Last  Updated:  28  Feb  97 
%  Written  By:  T.  H.  Newman 
%  Operating  Environment:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description:  Kalman  filter  for  tracking  phase  of  a 
%QPSK  signal 
% 

%  Based  on  Professor  Titus'  algorithim  which  follows  Prof.  Kirk’s 
%  NPS  lecture  notes "  Optimal  Estimation:  An  Introduction  to  the 
%  Theory  and  Applications"  1975  All  equation  numbers  are  referenced 
%  to  the  lecture  notes. 


function  y=kalman(z) 

global  Pkkml ;  %  Values  initialized  in  kal  init 

global  R; 

global  Q; 

global  I; 

global  Phi; 

global  xkkml; 

global  w; 


C=[1,0]; 

Cps=[0,l]; 

G=Pkkml*C'*inv(C*Pkkml*C'  +R); 
l)*C’*inv(C*P(k|k-l)*C'  +  R) 


%  Calculate  Gain  at  time  k 

%  G(k)  =  P(k|k- 


%  eqn4.3-102a 


Pkk=a-G*C)*Rdcml ;  %  Calculate  P(k|k) 

%  P(k|k)  =  [I  -G(k)*C(k)]*P(kik-l)  eqn  4.3-103a 

Pkkml=Phi*Pkk*Phi’  +  Q;  %P(klk-l) 

eO  =  z-C*xkkml;  xO=xldanl+G*eO; 

el  =  z-Cps*xkkml/w;  xl=Cps*(xkkml+G*el)/w; 
e2  =  z+C*xkkml;  x2=-(xkkml+G*e2); 

e3  =  z+Cps*xkkml/w;  x3=  -Cps*(xkkml+G*e3)/w; 

if  abs(eO)<=abs(el); 
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e=eO;  %  Smallest  error  is  with 

xkk  =  xkkml  +G*(e);  %  old  phase 

xkkps  =  xkk; 
else; 

e=el ;  %  Smallest  error  is  +90 

xkk  =  xkkml  +G*(e);  %  is 
xklq>s  =  Cps*xkk/w; 
end 

if  abs(e2)  <=  abs(e) 
e=e2; 

xkk  =  xkkml  +  G*(e); 

xkkps  =  -xkk;  %  Smallest  error  is  -90 

end; 

if  abs(e3)  <=  abs(e) 
e=e3; 

xkk  =  xkkml  +  G*e;  %  Smallest  error  is  180 

xkkps  =  -Cps*xkk/w; 
end; 


xkkml=Phi*xkk;  %  eqn  4.3-89a 

yl=  C*xkk;  %  report  the  non  shilled  version 

y=[yl(l,l);  xkkps(l);  x0(l);  xl(l);  x2(l);  x3(l)]; 


%  End  of  File  kalman.m 


H.  KALMAN_U.M 


%  File  Name:  kalman_u.m 
%  Last  Updated:  26  May  97 
%  Written  B>’:  T.  H.  Newman 
%  Operating  Environment:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description:  Kalman  filter  for  tracking  phase  of  a 
%  unbalanced  QPSK  signal 
% 

%  Based  on  Professor  Titus'  algorithim  which  follows  Prof  Kirk's 
%  NPS  lecture  notes  "  Optimal  Estimation:  An  Introduction  to  the 
%  Theory  and  Applications"  1975  All  equation  numbers  are  referenced 
%  to  the  lecture  notes. 
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function  y=kalman(z) 


global  Pkkm  1 ;  %  Values  initialized  in  kal_init 

global  R; 

global  Q; 

global  1; 

global  Phi; 

global  xkknil; 

global  w; 

global  Al; 

global  A2; 


C=[1,0]; 
Cps  =  [0,1]; 


G=Pkkml*C’*inv(C*Pkkml*C'  +R); 
l)*C'*inv(C*P(k|k-l)*C'  +  R) 


%  Calculate  Gain  at  time  k 

%  G(k)  =  P(k|k- 

%  eqn  4.3-102a 


Pkk=(I-G*C)*Pkknil;  %  Calculate  P(k|k) 

%  P(k|k)  =  [I  -G(k)*C(k)]*P(k|k-l)  eqn  4.3-103a 

Pkkml=Phi*Pkk*Plii'  +  Q;  %  P(k|k-1) 

F=(l/A2);  %  Factor  for  unbalanced  quadature  phase  term 

cos_wt=C*xkkml ; 
sin_wt=-Cps*xkkm  1/w; 


eO  =  z-(cos_'wt  +  sqrt(2)*sin_wt); 
el  =  z-(cos_nt  -  sqrt(2)*sin_wt); 
e2  =  z+(cos_wt  -  sqrt(2)*sin_wt); 
e3  =  z+(cos_wt  +  sqrt(2)*sin_wt); 

if  abs(eO)<=abs(el); 

e=e0;  %  Smallest  error  is  with 

xkk  =  xkkml  +G*(e);  %  old  phase 

xkkps  =  (cos_wt+sqrt(2)*sin_wt); 
else; 

e=el ;  %  Smallest  error  is  +90 

xkk  =  xkkml  +G*(e);  %  is 
xkkps  =  ( cos_wt  -  sqrt(2)*sin_wt ); 
end 


if  abs(e2)  <=  abs(e) 
e=e2; 

xkk  =  xkkml  +  G*(e); 

xkkps  =  (cos_wt  -  sqrt(2)*sin_wt);  %  Smallest  error  is  -90 

end; 

if  abs(e3)  <=  abs(e) 
e=e3; 
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xkk  =  xkkml  +  G*e;  %  Smallest  error  is  180 

xkkps  =  (cos_wt  +  sqrt(2)*sin_w't); 
end; 


xkkml=Phi*xkk;  %  eqn  4. 3-89a 

yl=  C*xkk;  %  report  the  non  shifted  version 

y=[yl(l,l);  xkkps(l)]; 

%  End  of  File  kalman_u.m 


L  KALMAN2.M 


%  File  Name:  kalman2.m 
%  Last  Updated:  3  May  97 
%  Written  By:  T.  H.  Newman 
%  Operating  Environment:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description:  Kalman  filter  for  tracking  phase  of  a 
%  BPSK  signal. 

% 

%  Based  on  Professor  Titus'  algoritfaim  which  follows  Prof  Kirk's 
%  NPS  lecture  notes "  Optimal  Estimation:  An  Introduction  to  the 
%  Theory  and  ./plications"  1975  All  equation  numbers  are  referenced 
%  to  the  lecture  notes. 

function  y=kalman(z) 

global  Pkkml ;  %  Values  initialized  in  kal  init 

global  R; 

global  Q; 

global  I; 

global  Phi; 

global  xkkml; 

global  w; 


C=[L0]; 

Cps  =  [0,l]; 

G=Pkkml*C'*inv(C*Rkkml*C'  +R); 
l)*C’*inv(C*P(k|k-l)*C'  +  R) 


%  Calculate  Gain  at  time  k 

%  G(k)  =  P(k|k- 

%  eqn  4.3-102a 
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Pkk=(I-G*C)*Pkkm  1 ;  %  Calculate  P(k|k) 

% 


P(k|k)  =  [I  -G(k)*C(k)]*P(klk-l)  eqn  4.3-103a 


Pkkml=Phi*Pkk*Phi'  +  Q;  %  P(k|k-1) 

eO  =  z-C*xkkml;  xO=xkkinl+G*(eO); 
e2  =  z+C*xkkml;  x2=-(xkknil+G*e2); 


if  abs(eO)<=  abs(e2); 

e=eO;  %  Smallest  error  is  with 

xkk  =  xkkml  +G*(e);  %  old  phase 

xkkps  =  xkk; 
else 
e=e2; 

xldc  =  xkkml  +  G*(e); 

xkkps  =  -xkk;  %  Smallest  error  is  -90 

end; 


xkkml=Phi*xkk;  %eqn4,3-89a 

yl=  C*xkk;  %  report  the  non  shifted  version 

y=[}'l(l,l);  xkkps(l)]; 


%  End  of  File  kalman2.m 


J.  MOD2.M 


%  File  Name:  mod2.m 
%  Last  Updated:  18  Oct  96 
%  Written  B>':  T.  H.  Newman 
%  Operating  Environment:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description:  Returns  the  modulo  2  value  of  input 
%  if  input  is  multiple  of  2  returns  0 
%  else  returns  1 
% 

% 


function  x  =  mod2(y) 


if  (y/2)*10  ==  round(y/2)*10 
x  =  0; 
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else 

x=l; 

end 


%  End  of  File  niod2.m 


K.  MPATH.M 


%  File  Name;  mpath.m 
%  Last  Updated;  28  Apr  97 
%  Written  By;  T.  H.  Newman 
%  Operating  Environment;  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description;  This  function  simulates  the  multipath  environment 

function  sm  =  mpath(s) 

global  delay 
global  alpha 
global  tau 

len=length(delay); 

if  ( len  >  1 ) 

delay  (2;len)  =  delay(l;len-l);  %  Shift  delay  vector 
delay(l)  =  s; 

sm  =  sum(delay(tau).*alpha); 

else 

delay(l)=s; 

sm=s; 

end; 


%  End  of  File  mpath.m 


L.  MPATH_AP.M 


%  File  Name;  mpath  ap.m 
%  Last  Updated;  4Jun97 
%  Written  Bj':  T.  H.  Newman 
%  Operating  Environment:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 
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%  Description:  Kalman  filter  for  multipath  with  a  priori 
%  information. 

clear;  %  clear  old  variables 


w=2;  %  =  2*pi*f 

Tf^l5;  %  finish  time sec 
dt=l; 

noise_fac=0;  %  change  scale  on  noise 

%multipath  a  priori  information 

alpha  =  .5;  %  delayed  signal  amplitude 

tau  =  10;  %  time  of  delay 

Amp  =  1 ;  %  amplitude  of  line  of  sight  signal 


A=[0  1  0  0;  -w^2  0  0  0;0  0  0  1;0  0  -w^2  0]; 
B=[0w^2  0  0]'; 

C=[l  0  1  0]; 

[Phi,Del]=c2d(A,B,dt); 

I=eye(4); 

Pkkml=le6*I; 

R=l; 

0=0. 1*1; 
kmax=Tfidt+l  ; 


switchl  =  3 ;  %  time  of  phase  changes 

S\\itch2  =  9; 

time=0; 

for  (i=l  :kmax)  %  generate  line  of  sight  signal 

t(i)=time; 
if  t(i)  <=  switchl 

los(i)  =  cos(w*t(i));  %  line  of  sight  direct  path  signal 

elseif  t(i)  >=  switchl  &  t(i)  <=  switch2 
los(i)  =  sin(w*t(i)); 
else 

los(i)  =  cos(w*t(i)); 
end; 

time=time+dt; 
end; 
tl=t; 
time  =  t; 
clear  t; 

for  i=  1  :length(los)“tau  %  generate  delayed  signal 
x(i)  =  los(i+tau);  %  Line  of  sight  signal  x(t) 

delay(i)  =  los(i)*alpha;  %  delayed  signal  alpha*x(t-tau) 
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z(i)  =  x(i)  +  delaj'(i);  %  z  =  x(t)  +  alpha*x(t-tau) 
t(i)  =  time(i+tau);; 
end; 

z  =  z+noise_fac*randn(size(z));  %  add  in  noise 

wt  =  w*time(tau+l); 
wtmtau  =  w*tinie(l); 
xkkml=[  cos(wt) 

-w*sin(wt) 

alpha*cos(wtmtau) 

-w*alpha*sin(wtnitau)] ; 


for  (i=l:length(z)) 

G=Hdanl*C'*inv(C*Hdanl*C'  +R); 

Pkk=(I-G*C)*Pkkml; 

Pkknil=Phi*Pkk*Phi’  +  Q; 
eO=z(l,i)-[l  0  1  0]*xkkml(:,i); 
el=z(l,i)-[0  -1/w  1  0]*xldanl(:,i); 
e2=z(l,i)-[0  -1/w  0  -l/w]*xkkml(:,i); 
e3=z(l,i)-[l  0  0  -l/w]*xkkml(:,i); 

%e_vect(:,i)  =  [  eO,  el,  e2,  e3  ]'; 

if  abs(eO)<=abs(el)  &  abs(eO)  <=  abs(e2)  %  abs(e0)<=abs(e3);  %t(i)<=sl; 

%if  (t(i)  <=  switchl) 
e^eO; 

xkk(:,i)  =  xkkml(:,i)  +G*(e); 
jddq5s(l,i)  =  xkk(l,i); 
xkl^(3,i)  =  xkk;(3,i); 

elseif  abs(el)  <=  abs(eO)  &  abs(el)  <=  abs(e2)  &  abs(el)<=abs(e3)  %  si  <  t  <  sl+phi 
%elseif  (switchl  <=  t(i)  &  t(i)  <=  switchl+tau*dt) 
e=el; 

xkk(;,i)=xkkml(:,i)+G*(e); 
xkkps(l,i)  =  -xkk(2,i)/w; 
xkkps(3,i)  =  xkk(3,i); 

elseif  abs(e2)<=abs(e0)  &  abs(e2)<=abs(el)  &  abs(e2)<=abs(e3)  %  sl+tau  <  t  <  s2 
%elseif  (switchl+tau*dt  <=  t(i)  &  t(i)  <=?  switch2) 
e=e2; 

xkk( :  ,i)=xkkml  ( :  ,i)+G*e; 
xldq>s(l,i)  =  -xkk(2,i)/w; 
xkkps(3,i)  =  -xkk(4,i)/w; 

elseif  abs(e3)<=abs(eO)  &  abs(e3)<=abs(el)  &  abs(e3)<=abs(e2) 

%elseif  (switch2  <=  t(i)  &  t(i)  <=  switch2+tau*dt) 

%  s2  <  t  <  s2+phi 

e=e3; 

xldc(:,i)=xkknil  ( :  ,i)+G*e; 
xklq3S(l,i)=xkk(l,i); 
xldq5s(3  ,i)=-xkk(4,i)/w; 
else 

e=eO; 
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xkk(;,i)  =  xkkinl(;,i)  +G*(e); 
xkkps(l,i)  =  xkk(l,i); 
xkkps(3,i)  =  xkk(3,i); 

end; 

xkkml(;,i+l)=Phi*xkk(:4); 

end; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%o^%o/o%%%%%%%%o^o^®^o^%o^% 

%%%%%%% 

figure; 

mpath=delay; 

los=x; 

figure; 

subplot(3,l,l); 

plotaios,V); 

ylabel('los'); 

subplot(3,lj2); 

plot(t,mpath,'b'); 

ylabel('mpath'); 

subplot(3,l,3); 

plot(t,z,'g'); 

ylabel('z'); 

figure; 

subplot(2,l,l) 

plot(t,los,'r',t,xkkps(  1 ,  OjIv.  '); 
title('los  &  xldq)s(l,:)  (estimate  = 

subplot(2,l,2); 

plot(t,mpath,'r',t,xklq)s(3, ;  ),'b-. '); 
title('mpath  &  xkkps(3,:)'); 


%  End  of  File  mpath  ap.m 
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M.  MY  GLB.M 


%  File  Name:  my_glb.m 
%  Last  Updated:  16  Feb  97 
%  Written  By:  T.H.  Newman 
%  Operating  Environment:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description:  Global  declaration  for  gps  simulation  and 
%  kalman  filter. 


global  C;  %  Values  initialized  in  kal_init.m 

global  Pkkml; 
global  R; 
global  Q; 
global  I; 
global  Phi; 
global  xkkml ; 

global  w; 

fl  =  l; 

To=l; 
w=2*pi*fl; 
phase_init  =  pi/6; 
num_save  =  100000; 


%  NOTE:  These  are  the  same  as  in  init.m 
%  and  are  placed  here  for  use  in  non-gps  model. 

%  Set  frequency  for  simulation 
%  Phase  offset  for  sine  wave  in  BPSK  modulator 
%  Number  of  points  to  save  to  workspace  for  plotting  later 


%  End  of  File  my_glb.m 


N,  SHIFT  REG.M 


%  File  Name:  shift_reg.m 
%  Last  Updated:  21  Oct  96 
%  Written  By:  T,  H.  Newman 
%  Operating  Environment:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description:  This  function  generates  the  shift  register 
%  linear  code  generators  where  R  =  the  input  values 
%  size  =  register  size.  The  register  is  shifted  right 
%  and  returned  in  a  new  vector 

function  r  =  shift_reg(R,  size) 

r  =  ones(l,size) 


%  End  of  File  shrft  reg.m 
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o 


SIG  GEN.M 


%  File  Name;  sig_gen.m 
%  Last  Updated:  1 1  Nov  96 
%  Written  By:  T.  H.  Newman 
%  Operating  Enviromnent:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description:  Returns  either  1  P(Y)  code  xor  with  data 
%  2  P(Y)  code 

%  or  3  C/A  code  xor  with  data 

%  the  number  is  the  sig_num  which  determines  which  signal 
%  is  returned 
% 


function  [sl,s2,s3]  =  sig  genfcode  length) 


sv=l; 

xlepoch  =  0; 
xlA_epoch  =  0; 
xlBepoch  =  0; 
x2A_epoch  =  0; 
x2B_epoch  =  0; 

%code_length  =  10; 

r_gl  =  ones(  size(l :  10) ) ;  %  initialize  register  values 

r_g2  =  ones(  size(l :  10) ); 

r_xla=[000  100  100  100]; 

r_xlb=[00  10  10  10  10  10]; 

r_x2a=[10  100  100  100  1]; 

r_x2b=[00  10  10  10  10  10]; 

shift_register  =  zeros(size(l:37)); 


t  =  0; 


dataclock  =  20; 

for  i  =  1  :code_length 
if  data_clock  ==  20 

data  =  1 ;  %  set  to  all  ones  to  ignore 

data_clock  =  0; 
else 


data_clock  =  data_clock  +  1 ; 

end; 

for  j  =  1 ;  10  %  gen  C/A  code  at  1 .023  MHz 

[r_gl,  g_l]  =  gl(r_gl,  xl_epoch); 
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[r_g2,  g_2]  =  g2(r_g2,  xl_epoch,  sv); 
g  =  xor(g_l,g_2); 


for  k=  1 ;  10  %  gen  P/A  code  at  10.23  MHz 

t=t+l;  %  or  10  *  C/A  code  generation  rate 

[r  xla,  A]  =  xla(r_xla,  xlA  epoch); 

[r_xlb,  B]  =  xlb(r_xlb,  xlB_epoch); 

[r_x2a,  C]  =  x2a(r_x2a,  x2A_epoch); 

[r_x2b,  D]  =  x2b(r_x2b,  x2B_epoch); 

xl=xor(A,B); 

x2=xor(C,D); 

shifl_register(2:37)  =  shift_register(l;36); 
shift_register(l)  =  x2; 
p=xor(xl,  shift_register(sv)); 

sl(t)  =  xor(p,data); 
s2(t)  =  p; 

s3(t)  =  xor(g,  data); 
end;  %  for  k; 

end;  %  for  j 

end;  %  for  i 


%  End  of  File  sig  gen  m 


P.  TEST_CA,M 


%  File  Name:  test_CA.m 
%  Last  Updated:  3 1  Mar  97 
%  Written  B>':  T.  H.  Newman 
%  Operating  Environment:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description:  Script  file  to  test  C/A  code  generation.  Output 
%  should  match  First  10  C/A-Chips  (Octal)  of  GPS-lCD-2000  table 

for  sv=l:37 

r_gl  =  roimd(  rand(  size(l :  10) ) );  %  initialize  to  some  random  state 
r_g2  =  round(  rand(  size(l:10) ) ); 
xl  =  1; 

i=i; 

[r _gl,  g_l(i)]  =  gl(r_gl,  xl);  %  Initialize  at  epoch 

[r_g2,  g_2(i)]  =  g2(r_g2,  xl,  sv); 
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%  Simulate  1st  10  chips 


xl=0; 

fori=2:10 

[r_gl,  g_l(i)]  =  gl(r_gl,  xl); 

[r_g2,  g_2(i)]  =  g2(r_g2,  xl,  sv); 
end; 

g  =  xor(g_l,g_2);  % 

digl  =  g(l); 

(lig2  =  g(2)*4  +  g(3)*2  +  g(4);  %  convert  to  octal  number 

dig3  =  g(5)*4  +  g(6)*2  +  g(7); 

dig4  =  g(8)*4  +  g(9)*2  +  g(10); 

result(sv)  =  digl*10''3  +  dig2*10^2  +  dig3*10  +  dig4; 

end 

result  %  display  results 

%  End  of  File  test_CA.m 


Q.  TEST_PJ\I 


%  File  Name:  test  P.m 
%  Last  Updated:  1  April  97 
%  Written  By:  T.  H.  Newman 
%  Operating  Environment:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description:  Script  file  to  test  P  code  generation.  Output 
%  should  match  First  12  P-Chips  (Octal)  of  GPS-ICD-2000  table 
%  Note  that  only  the  first  13  sv's  are  calculated.  SV  9  -  37  are 
%  identical 

for  sv=  1:13 

r_xla  =  round(  rand(  size(l:  12) ) );  %  initialize  to  some  random  values 
rxlb  =  round(  rand(  size(l :  12)  )  ); 
r_x2a  =  round(  rand(  size(l:  12) ) ); 
r_x2b  =  round(  rand(  size(l:12) ) ); 

epoch  =  1; 
i=l; 

[r  xla,  A(i)]  =  xla(r_xla,  epoch); 

[r_xlb,  B(i)]  =  xlb(r_xlb,  epoch); 

[r_x2a,  C(i)]  =  x2a(r_x2a,  epoch); 

[r_x2b,  D(i)]  =  x2b(r_x2b,  epoch); 

epoch=0; 

fori=2:12 

[r_xla,  A(i)]  =  xla(r_xla,  epoch); 
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[r_xlb,  B(i)]  =  xlb(r_xlb,  epoch); 

[r_x2a,  C(i)]  =  x2a(r_x2a,  epoch); 

[r_x2b,  D(i)]  =  x2b(r_x2b,  epoch); 
end; 

xl=xor(A3); 

x2=xor(C,D); 

delay=sv; 

xl=[xl,  zeros(l, delay)]; 
x2=[ones(l, delay),  x2  ]; 
p=xor(xl,  x2); 

digl  =  p(l)*4  +  p(2)*2  +  p(3); 
dig2  =  p(4)*4  +  p(5)*2  +  p(6); 
dig3=p(7)*4  +  p{8)*2+p(9); 
dig4  =  p(10)*4  +  p(l  1)*2  +  p(12); 
result(sv)  =  digl*10''3  +  di^*10''2  +  dig3*10  +  dig4; 
end;  %  for  sv 
result 

%  End  of  File  test  P.m 


R.  UI_INIT.M 


%  File  Name:  ui_init.m 
%  Last  Updated:  4  May  97 
%  Written  By:  T.  H.  Newman 
%  Operating  Environment:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description:  Called  by  ui_vars.m  to  initialize  parameters 
%  after  Ae  Init  button  is  pushed. 

% 


delay=zeros(size(l:max(tau) )); 
if  SNR  >0 

Noise _gain  =  sqrt((Amp^2)/(  10''(SNR/10) )); 
else 

Noise _gain  =  1; 
end; 

if  SNR  >=100 
Noise_gain  =  0; 
end; 

w=2*pi*f; 


disp(sprintf('  f=%2.2f  Initial  Phase  =  %2.2f ,  f,  p  degree)) 
disp(sprintf('  Amp=%2.2f,  Amp)) 
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disp{sprintf('  dt=%2.2f,  stop  =  %2.2f,  Tb=%2.2f ,  dt,  stop,  Tb)) 
disp(sprintf('  SNR=%2.2f  dB,  Noise_gain=%2.2f  ,SNR,  Noise_gain)) 
disp('tau='),  disp(tau) 
disp('alpha='),disp(alpha) 


%  End  of  File  ui  initm 


S,  UI_VARS.M 


%  File  Name:  ui_vars.m 
%  Last  Updated:  28  Apr  97 
%  Written  By:  T.  H.  Newman 
%  Operating  Environment:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description:  An  Editable  Text  windows  to  control 
%  initialization  values  for  kal_init 
%  Calls:  my_global.m  to  initialize  values 
%  kal_init.m  to  initialize  kalman  filter  variables 

%  init.m  to  generate  GPS  code  sample 

global  MainFig 

%  below  are  the  initial  values  that  will  be  in  the  text  window 

global  C;  %  Values  initialized  in  kal_init.m 

global  Pkkml; 

global  R; 

global  Q; 

global  I; 

global  Phi; 

global  xkkml; 

global  Tb; 

global  Noise_^n; 

global  dt; 

global  Amp;  %  Amplitude  for  BPSK 

global  tau; 

global  alpha; 

global  delay; 

global  w; 

global  clock; 

my_glb; 

%  This  initilizes  w,  phase  init,  num  save,  Tb,  Noise  gain,  Al,  A2,  dt 
%  Retreive  initial  values  from  the  strings 
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start  =  0;  %start  time 

stop  =  10;  %stop  time 

f  =  w/(2*pi);  %signal  frequency 

p_degree  =  phase_init*  180/pi;  %phase  in  degrees 

SNR  =100; 

f_string  =  sprintfC  %2.2f,f);; 
p_degree_string  =  sprintfC  %2.2f  ,p_degree); 
Amp_string  =  sprintfC  %2.2f  ,Amp); 

dt_string  =  sprintfC  %2.2f  ,dt);; 

start_strmg  =  sprintfC  %2,2f, start); 

stopstring  =  sprintfC  %2.2f  ,stop); 

Tb_string  =  sprintfC  %2.2f  ,Tb); 

SNR_string  =  sprintfC  %2.2f,  SNR); 

tau_string  =  '  [1  ]'; 

alpha_string  = '  [1  ]'; 

tau=eval(tau_string); 
alpha  =  eval(alpha_string); 

%define  the  Main  Figure  Window. 

num_entries  =  9; 

sep=45; 

ymax  =  sep*(num_entries)+10; 

x_pos=25;  y_pos=35;  dx=125;  dy  =  ymax; 

MainFig  =  figure  (Position',  [xjx)s  y_pos  dx  dy],... 
'Color',  'white','NumberTitle','off,'Name',. . . 
'Simulation  Variables',... 

'MenuBar','none'); 


o/„ - ; - % 

%  Define  UIMENU  for  pull  down  menu  options 

e=uimenuClaber,'liiit', 'callback', . . . 

'kal_init,ui_init; '); 

f=uimenuClaber,  Plot'); 

oimenu(f,  'laber,Phase  Lock  Loop','callback','plot_k'); 

uimenu(f,1aber,'Delay','callback','plot_d'); 

uimenu(f,'laber,'QPSK','callback','plot_q2'); 

uimenu(f,'label',PPSK','callback','plotbpsk'); 

uimenu(f,1aber,'Error','callback','plot_err'); 

uimenu(f,'laber,'Variables','callback','plot_k2'); 


g=uimenuClabel',  'Sim'); 
gl=uimenu(g,  'laber,'BPSK'); 
uimenu(gl,  'label',  'BPSK  Simulation','caIlback','kal_sim') 
uimenu(gl,  'label',  'PLL  Simulation','caUback','pir); 
g2=uimenu(g,  'label',  'QPSK'); 


81 


uimenu(g2,  'label'/QPSK  Simulation',  'callback','kalq_sim'); 
uimenu(g2,  'label',  'PLL  Simulation','callback','pll_q'); 
g3=uimenu(g,  'label'.  Unbalanced  QPSK'); 
uimenu(g3,  'label', 'Unbalanced  QPSK  Simulation',  'callback', 'kalu_sim'); 
uimenu(g3,  'label',  PLL  Simulation','callback','pll_u'); 
g4=uimenu(g,'laber,'Multipath','callback','mpathsim'); 


n=l;  %uimenu  number 

jinax  =  dy; 

xpos=10; 

dx=  100; 

dy=20; 

space  =  22; 

label_fg_color  =  'Black'; 
labelJbg_color  =  'White'; 
boxfgcolor  =  'Black'; 
box_bg_color  =  'Yellow'; 

f_label  =  uicontrol(MainFig,  'Style',text',  'Pos',[xpos  ymax-n*sep+space  dx  ^],... 
'String',  'Signal  freq(Hz)',PoregroundColor',label_fg_color,... 
'BackgroimdColor',label_bg_color); 


f_box  =  uicontrol(MainFig,’Style','edit','String',  f_strmg,... 
'Pos',[xpos  ymax-n*sep  dx  dy],... 

'Foregroundcolor', box_fg_color,... 
'BackgroundColor',box_bg_color,... 

'Callback','f  =  eval(  get(f_box,"String")  );  w=2*pi*f;'); 


n=n+l; 

phase_label  =  mcontrol(MainFig,  'Style', 'text',  'Pos',[xpos  ymax-n*sep+space  dx  <fy],... 
'String',  'Phase  Offset  °',PoregroundColor',label_fg_color,... 
'BackgroundColor',label_bg_color); 

phase_box  =  uicontrol(MainFig, 'Style', 'edit','String',  p_degree_string,... 

'Pos',[xpos  ymax-n*sep  dx 
Toregroimdcolor',box_fg_color, . . . 

'BackgroundColor',box_bg_color, . . . 

'Callback','p_degree  =  eval(  get(phase_box,"String")  );phase_mit=p_degree*pi/180'); 
n=n+l; 

amp  label  =  uicontrol(MainFig,  'Style','text',  Pos',[xpos  ymax-n*sep+space  dx  <ty],... 
'String',  'Amplitude','ForegroundColor',label_fg_color,... 
'BackgroundColor',label_bg_color); 

amp_box  =  uicontrol(MainFig,'Style','edit','String',  Amp_string,.., 

'Pos',[xpos  ymax-n*sep  dx  dy],... 
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'Foregroundcolor', boxfgcolor,. . . 

'BackgroundColor',box_bg_color, . . . 

'Callback','Al  =  eval(  get(amp_box, "String") );'); 

n=n+l; 

sanip_label  =  uicontrol(MainFig,  'Style','text',  'Pos',[xpos  ymax-n*sep+i5pace  dx  dy], 
'String',  'Ts  (sec)','ForegroundColor',label_fg_color,... 
'BackgroundColor',label_bg_color); 

samp_box  =  uicontrol(MainFig,'Style','edit',’String',  dt_string,... 

'Pos',[xpos  ymax-n*sep  dx 
Toregroundcolor',box_fg_color,. . . 

'BackgroundColor',box_bg_color, . . . 

'Callback','dt  =  eval(  get(samp_box, "String") );'); 

n=n+l; 

stop_label  =  uicontrol(MainFig,  'Style','text',  'Pos',[xpos  yniax-n*sep+space  dx  dy],.. 
'String',  'Stop  time  (sec)','ForegroundColor',label_fg_color,. . . 
'BackgroundColor',label_bg;_color); 

stop_box  =  uicontrol(MainFig,'Style','edit','String',  stop_string,... 

'Pos',[xpos  ymax-n*sep  dx  dy'],... 

'Foregroundcolor', box_fg_color,. . . 

'BackgroundColor',box_bg_color, . . . 

'Callback’,'stop  =  eval  ( get(stop_box,"String") );'); 

n=n+l; 

Tb_label  =  uicontrol(MainFig,  'St5'le','text',  'Pos',[xpos  ymax-n*sep+space  dx  ^],... 
'String',  'Tb  (sec)','ForegroundColor',label_fg_color,... 
'BackgroundColor',label_b^color); 

Tb_box  =  uicontrol(MainFig,'Style','edit','String',  Tb_string,... 

'Pos’,[xpos  yniax-n*sep  dx  <fy],... 

'Foregroundcolor',box_fg_color, . . . 

'BackgroundColor',box_bg_color, . . . 

'Callback','Tb  =  eval  ( get(Tb_box,"String") );'); 

n=n+l; 

SNR_label  =  uicontrol(MainFig,  'Style', 'text',  'Pos',[xpos  ymax-n*sep+space  dx  <fy],.. 
'String',  'SNR  dB','ForegroundColor',label_fg_color,... 
'BackgroundColor',label_bg_color); 

SNR_box  =  uicontrol(MainFig,'Style','edit','String',  SNR_string,... 

'Pos',[xpos  ymax-n*sep  dx  ^'],... 

Toregroundcolor',box_fg_color,. . . 

'BackgroundColor',box_bg_color, . . . 

'Callback','SNR  =  eval(  get(SNR_box, "String") ); '); 
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n=n+l; 


tau  label  =  uicontrol(MainFig,  'Style','text',  'Pos',[jq)os  ymax-n*sepfspace  dx 
’String',  'tau','ForegroundColor’,label_fg_coIor,. . . 
'BackgroxmdColor',label_bg_color); 

tau_box  =  uicontrol(MainFig,  'Style', 'edit', 'String',  tau_string,... 

'Pos',[3^s  ymax-n*sep  dx  dy],... 

'Foregroundcolor',box_fg_color,. . . 

'BackgroundColor',box_bg_color, . . . 

'Callback', 'tau  =  eval(  get(tau_box,"String") ); '); 

n=n+l; 


alphajabel  =  uicontrol(MainFig,  'Style', 'text',  'Pos’,[jq)os  ymax-n*sepfspace  dx  cfy],... 
'String',  'alpha', 'ForegroundColor',label_fg_color,... 
'BackgroundColor’,label_bg_color); 

alpha_box  =  uicontrol(MainFig,'Style','edit', 'String',  alpha_string,... 

'Pos',[xpos  yniax-n*sep  dx  dy],... 

'Foregroundcolor',box_fg_color, . . . 

’BackgroundColor',box_bg_color, . . . 

'Callback', 'alpha  =  eval(  get(alpha_box,"Striing") ); '); 

n=n+l; 

O/o - o/^ 

%  Define  the  UICONTROL  button  for  callbacks. 

%  bl=mcontrol('style',’pushbutton','string','Initiali2e', . . . 

%  'position',  [30  ymax-n*sep  65  25],'Callback',... 

%  'kalinit;  delay=zeros(size(l;max(tau) ));'); 


%  End  of  File  ui  vars.m 


T.  XIAJM 


%  File  Name:  xla  m 
%  Last  Updated;  18  Oct  96 
%  Written  By:  T.  H.  Newman 
%  Operating  Environment:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 
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%  Description:  This  function  generates  the  xla(t)  signal 
%  for  GPS  P  Code  generator 

function  [new_R,  A]  =  xla(R,  xla_epoch) 

if  xla_epoch==l 

R  =  [0  0  0  1  0  0  1  0  0  1  0  0];  %  Reset  at  XI A  epoch 
end; 

new_R  =  gen_poly([l, 6,8,1 1,12],  R); 

A=R(12); 

%  End  of  File  xla.m 


U.  X1B.M 


%  File  Name:  xlb.m 
%  Last  Updated:  18  Oct  96 
%  Written  By:  T.  H.  Newman 
%  Operating  Environment:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description:  This  function  generates  the  xlb(t)  signal 
%  for  GPS  P  Code  generator 

function  [new_R,  B]  =  xlb{R,  xlb_epoch) 


if  xlb_epoch==l 

R=  [00  1 0  1  0  1  0  1  0  1  0];  %ResetatXlB epoch 
end; 

new_R  =  gen_poly([l  12  5  8  9  10  11  12],  R); 
B=R(12); 

%  End  of  File  xlb.m 


V.  X2AJVI 


%  File  Name:  x2a.m 
%  Last  Updated:  18  Oct  96 
%  Written  By:  T.  H.  Newman 
%  Operating  Environment:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description:  This  function  generates  the  x2a(t)  signal 
%  for  GPS  P/A  Code  generator 
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fimction  [new  R,  C]  =  x2a(R,  x2a_epoch) 


if  x2a_epoch==l 

R=  [1  0  1  GO  1  0  0  1  0  0  1];  %  Reset  at  X2a  epoch 
end; 

new_R  =  gen_poty([l  1  3  4  5  7  8  9  10  11  12],  R); 
C=R(12); 

%  End  of  File  x2a.m 


W.  X2B.M 


%  File  Name:  x2b.m 
%  Last  Updated:  18  Oct  96 
%  Written  By:  T.  H.  Newman 
%  Operating  Environment:  Windows  95 
%  Matlab  Version  4.2c.  1 
% 

%  Description:  This  fimction  generates  the  x2b(t)  signal 
%  for  GPS  P  Code  generator 

fimction  [new_R,  D]  =  x2b(R,  x2b_epoch) 


if  x2b_epoch==l 

R=  [00  1  0  1  0  1  0  1  0  1  0];  %  Reset  at  X2B  epoch 
end; 

new_R  =  genjx)ly([l  2  3  4  8  9  12],  R); 

D=R(12); 

%  End  of  File  x2b.m 
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APPENDIX  B  SIMULINK  DIAGRAMS 


A.  BPSK.M 


Outport 


Figure  35  BPSK.M 


-►t  BPSK  I 
To  Workspace 
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BPSK  P90.M 


Outport 


Figure  36 


BPSK_P90.M 
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BPSK  U.M 


Outtaort 


Figures?  BPSKJJM 


D.  KAL_SIM.M 


Number 


Kalman  Filter  as  Phase  Lock  Loop  Simulation 


Figure  38  KAL  J1M.M 


E.  KALQ_SIM.M 

O - W  t  I 

Clock  time 


Random 

Number! 


Figure  39  KALQ  SIMM 


F.  KALU_SIM.M 

© - 

Clock  time 


Numberl 

Kalman  filter  Phase  Lock  Loop  Simulation 
w/  Unbalanced  QPSK  signal 
(4  possible  phase  changes) 


Figure  40  KALUJIMM 


G.  LI  SIM 


© — 

Clock  tima 


I  sig3  |> 
C/A  xor  Data 


Note:  Run  initm  file  first  to  set  parameters 
and  generate  sig1,  sig2,  sig3 


Figure  41  L1_SIM.M 
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'  sigi  ' 


PY  xor  Data 


p!l 


S-function2 


sig3  - 

— ►  bpsk_p90  ■ 

C/A  xor  Data 

S-functionI 

Note:  Run  initm  file  first  to  set  parameters 
and  generate  sigi ,  sig2,  sig3 


Figure  42  L2_SIM.M 
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L 


MPATHSIM.M 


Figure  43  MPATHSIM.M 
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J. 


PLL.M 


Phase  Lock  Loop 
pll.m 

Figure  44  PLL.M 
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K.  PLL  2A.M 


Demux 


Figure  45  PLL2A.M 
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APPENDIX  C  EXAMPLE  SNR  CALCULATION 


dB2F(x) :  =  10^®  convert  dB  to  ratio  dB(x)  :=  lO-log(x)  convert  ratio  to  dB 

Minimum  Received  Power  C  =  - 160.0  ~  dB_W 

L1  with  C/A  Code 


Boltzman's  Constant 


k:=- 228.6  ~dBW/Hz-K 


Noise  Temperature 
Received  Noise 


Tn  :=28  ~  dB  K 

No  :=k+  Tn 


Carrier  to  Noise  Density  Ratio  CNR C  -  No 

CNR=40.6  ~dbHz 


Chip  Data  Rate 

R:=  1.023 10^ 

~  chips/sec 

R_db  :=dB(R) 

R_db  =60.099 

~  dB  cps 

Information  Data  Rate 

II 

o 

~  bits/sec 

Band  Width  -  QPSK 

B:=0.6R 

B_dB:=dB(B) 

B_dB  =  57.88 

~  dB  bps 

Signal  to  Noise  Ratio 

Processing  Gain 
Signal  to  Noise  Ratio 


SNR  :=CNR- B_dB 
SNR=-17.28  ~  dB 

/  T?  \ 

PG:=dB  — )  PG=43.109  ~  dB 

\R_i/ 

SNR:=SNR-rPG  SNR= 25.829  ~dB 
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